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The formalism of Kundu etal. [J. Stat. Mech. (201 1) P03007], for computing the large deviations of heat flow 
in harmonic systems, is applied to the case of single Brownian particle in a harmonic trap and coupled to two heat 
baths at different temperatures. The large-T form of the moment generating function (e^^^) ~ g{X) exp[t fl{X)], 
of the total heat flow Q from one of the baths to the particle in a given time interval T, is studied and exact 
explicit expressions are obtained for both /J (A) and g(A). For a special case of the single particle problem that 
corresponds to the work done by an external stochastic force on a harmonic oscillator coupled to a thermal bath, 
the large-T form of the moment generating function is analyzed to obtain the exact large deviation function as 
well as the complete asymptotic forms of the probability density function of the work. 
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I. INTRODUCTION 

The emergence of the so-called fluctuation relations |[T}j7l 
has generated considerable theoretical |8- 19| and experimen- 
tal 120-33 1 interests in studying fluctuation of various stochas- 
tic quantities such as entropy production, heat flow, particle 
transfer, power injection, work done, etc. in a given time 
interval T, in nonequilibrium systems. In this context, one 
is usually interested in the large-T behavior of the probabil- 
ity distribution of the studied stochastic quantity, say Q. The 
probability distribution P{Q) is expected to have a large devi- 
ation form 1 34 1 P{Q) ^ exp[Th{Q/T)] for large T. However, in 
spite of the large interests, there are only a few systems where 
the large deviation function h{q) has been obtained exactly. 

An interesting example is a free Brownian particle coupled 
to two heat baths at different temperatures, initially introduced 
by Derrida and Brunei |35 1. hi this model, the quantity of in- 
terest is the heat flow from one of the baths to the particle 
in a given time interval T. It turns out that this model can 
be mapped to an exactly solved problem, namely, the quan- 
tum harmonic oscillator — Visco ifTTll has obtained the ex- 
act characteristic function (e^^^) for all T, and thereafter, the 
large deviation function h{q) by analyzing the large-T limit of 

Now in addition, if a harmonic trap is introduced around 
the Brownian particle, this seemingly simple model becomes 
quite non-trivial and any methods for solving the relevant 
Fokker-Planck equation [Eqs. (|7| and ([8])] is not known to 
us. Fortunately, one does not require the complete solution 
of the Fokker-Planck equation, as the large-T behavior is es- 
sentially determined by the largest eigenvalue and the coiTe- 
sponding eigenfunctions (left and right) of the Fokker-Planck 
operator [see Eq. (|9]l]. It is indeed remarkable that these func- 
tions could be found exactly (in terms of some integrals) for 
a harmonic chain whose two ends are coupled to heat baths at 
different temperatures |36|. The harmonically bound Brow- 
nian particle is a special case of the harmonic chain. In this 
paper we find the eigenvalue and the eigenfunctions explicitly 
for the harmonically bound particle, and using those we then 
find the characteristic function (e^^^) for large T. In fact, a 
special case of this model — that concerns the fluctuations of 
the work done by an external stochastic force on a harmonic 



oscillator coupled to a thermal bath — has been realized in a 
recent experiment ||33l . The goal of the this paper is to ana- 
lyze this problem in detail, as the methods should be useful for 
other similar problems. Some of the main results have been 
reported in |37|. 

This paper is organized as follows. In Sec.|ll]we discuss the 
problem of the harmonically bound Brownian particle coupled 



to two heat baths. Section III contains the special case of the 



harmonic oscillator driven by a random force. In Sec. Ill A 
we obtain the large deviation function of the work fluctua- 



tion and in Sec. Ill B we find the complete asymptotic form of 
the probability density function of the work. Section IV con- 



tains some remarks on the evaluation of stochastic integrals 
that define heat, work, etc., in numerical simulations or from 
experimental data. Finally we summarize in Sec. [V] Some of 
the details are presented in the appendices. In AppendixjAjwe 
outline the relevant results of Ref. [36 1 for the harmonic chain 



and Appendix A 1 contains some of the details of the bounded 
Brownian particle case discussed in Sec. |ll] In Appendix |B] 
we outline the method of uniform asymptotic expansion of a 
integral having a saddle-point near a pole, which is used to ob- 
tain the complete asymptotic form of the probability density 
function in Sec. lIIIBl 



II. BROWNIAN PARTICLE IN A HARMONIC 
POTENTIAL COUPLED TO TWO THERMOSTATS 



Consider a Brownian particle of mass m in a one- 
dimensional harmonic potential with spring constant k and 
coupled to two white noise Langevin thermal baths at two dif- 
ferent temperatures and respectively. The displacement 
x{t) from its mean position and the velocity v(f ) of the particle 
are described by the equations 



and mv — —yv — kx + r]{t), 



(1) 



where 7 = 7l + 7ff and T7(f) = I7l(0 + m{t) with rji and t]^ 
being the Gaussian white noises with mean zero and correla- 
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tors: 

{riL{s)riLit)) - 2dL5{s - 1) with = YtkBTL, (2) 

{riR{s)riR{t)) ^2dR5{s-t) with dR ^ YRkBTn, (3) 

and {riL{s)riRit))=0, (4) 

where kg being the Boltzmann constant. The quantity of in- 
terest is the total amount of heat flowing from the reservoir at 
temperature Ti to the particle in a time duration t: 



Q 



[r]L{t)-rLvit)]v{t)dt. 



(5) 



This system is a special case (A^ = 1) of the harmonic chain 
(of particles) connected at its two ends to reservoirs at dif- 
ferent temperatures, that has been studied recently ll36l and is 
outlined in Appendix [A] 

We consider the restricted characteristic function 



Z{X,x,v,T\xo,vo)^{e-^Q8[x-x{T)]5[v-v{T)]) 



(6) 



for fixed initial and final configurations, (jcq, vq) and {x,v) re- 
spectively. It satisfies the Fokker-Planck equation j 



Z(A,x,v,t|xo,vo) = ifAZ(A,x,v',T|xo,vo) 



(7) 



with the initial condition Z(A,ji:,v,0|xo,vo) — — xo)5(v — 
vq) and the Fokker-Planck operator is given by 



-dR 



k y+2ldi 

— xH V 

m m 



d 



d 

dx 



-,-+xij,+xd,y+^-±^. 



(8) 



The solution of the Fokker-Planck equation can be formally 
expressed in the eigenbases of the operator and the large 
T behavior is dominated by the term having the largest eigen- 
value. Thus, for large T, 



Z{X,x,v,x\xq,vq) --x(xo,vo,A)1'(jc,v,A)e 



T/j(l) 



(9) 



where ^(x,v,A) is the eigenfunction corresponding to the 
largest eigenvalue /i(A) and x{xq,V(),X) is the projection of 
the initial state onto the eigenstate corresponding to the eigen 
value /i(A). These functions are obtained in Appendix 
and we find that 



Al 



^ix,v,X) 



exp[-B+(A)£'(x,v)], 



2t. 

Yri{X)y/km 
2n(dL + dR) 

and x(xo,vo,A) =exp[-B_(A)£'(jco,vo)], 

where Ty = m/7 is the viscous relaxation time. 




l+4^A(Aj3-A) 



kx^ 



(10) 

(11) 
(12) 

(13) 
(14) 



andB±(A) is given by Eq. (A27i. Fogedby and Imparatohave 
recently shown Ii38i that /i(A) can also be obtained by the 
Derrida-Brunet method ||35]| . 

Using the explicit forms one can easily verify that 
^X^{x,v,X) = jU(A)^(x, v,A). Moreover, since from \A21\ 
we get B+(A)+B_(A) = 'Yri{X)/ {dL + dR), it immediately 
follows that 



/:/: 



X{x,v,XYi'{x,v^X)dxdv — 1, 



(15) 



which is demanded by the normalization. From the above 
expressions, we also find that ju(0) =0 and ;if(xo,vo,0) = 1. 
Since A = case of Eq. ^ gives the probability distribution of 
the phase-space variables and ji (0) is the largest eigenvalue, it 
follows from Eq. (j9| that ^(x, v, 0) is the steady-state distribu- 
tion of the phase-space. Therefore, averaging over the initial 
variables (xo,v'o) with respect to ^(xo,vo,0) and integrating 
over the final variables (x,v), we find the characteristic func- 
tion of the heat flow in the steady state as 



where 



Z{X,x) = {e-^Q)^g{X)e'^(^\ 
47] (A) 



'\ + n{Xf-[2XdUyf' 



(16) 



(17) 



Interestingly, both /x(A) and g{X) are independent of the 
spring constant k. However, while /i(A) is same for both k^Q 
and /t = cases (the latter was obtained in Ref. [ 17|), ^(A) for 
k^O differs from that for k — Q. More precisely, g{X)\j^^Q = 
[g(A)|i.=o]^. The k^Q limit of ^(A) is not same as the k — 
case. Therefore, although the large deviation function are 
same for both the cases, the precise asymptotic form of the 
probability density functions of Q are different. 

The leading behavior of the probability density function 
P{Q) ^ exp[T/i(2/T)] can be obtained by inverting Eq. ( 16 1 
using the saddle-point approximation. Usually the prefactor 
g{X) can be ignored in such calculation and the large deviation 
function h{q) is related to /i(A) by the Legendre transform 



h{q)=^{X*)+X*q, 



-H'{X*)=q. 



(18) 



However, if g{X) has any singularities in the region of the 
saddle-point integration, the functions h{q) and /i(A) are not 
simply related by the Legendre transform, and it it impor- 
tant to retain the prefactor ^(A) in the saddle-point calcula- 
tion Il8l|9l[l71, as we see in the next section where we consider 
a special case that corresponds to an experiment reported in 
Ref [ITl. 



III. HARMONIC OSCILLATORS DRIVEN BY AN 
EXTERNAL RANDOM FORCE 

Consider a harmonic oscillator coupled to a thermal bath 
and driven out of equilibrium by an external Gaussian random 
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force. The displacement x{t) of the harmonic oscillator from 
its mean position is described by the Langevin equation 



d X 



dx 



(19) 



where m is the mass, 7 is the viscous drag coefficient and k 
is the spring constant. The interaction with the thermal bath 
is modeled by a Gaussian white noise CrlO with zero-mean 
{Crit)) — 0- The externally applied force /o(f) is again a 
Gaussian random variable with (/o(r)) — 0' ™d and /o 
are uncorrected. Equation ( 19 1 is asymmetric in and /o — 



the fluctuation-dissipation theorem relates the thermal fluctu- 
ation to the viscous drag as {tl,j[s)tl,j[t)) = 2D5{s — t) where 
D — yksT with T being the temperature of the bath and kg 
being the Boltzmann constant, whereas the fluctuation of the 
external force {fo{s)fo{t)) = {5fo)^5{s — t) is independent of 
7. The quantity of interest is the work done by the external 
random force /o(f) on the harmonic oscillator in a time in- 
terval T, in the nonequilibrium steady state. This is given (in 
units of ksT) by 



1 



(20) 



with the initial condition (at T = 0) drawn from the steady 
state distribution. 

It is evident that this harmonic oscillator problem can be 
mapped to the problem of the Brownian particle discussed 
in the previous section, with the following set of transforma- 
tions: 



• /o(f) with Yl 

YlTl = di fixed and di^ - 



^■ 0, oa while keeping 
(5/o)V2. 



• Tjff(f) ^rit) with Tr -^T,Yr^ 7, and dR D. 

Under these transformations, we have Q {y/Dy^Wr. 
Therefore, using Eq. ( [T6| we now get 



(21) 



where and g(A) are given by Eq. ( 10 1 and Eq. ( 17 1 re- 
spectively with the transformation X (j/D) X. 

It is useful to express the relative strength of the external 
force with respect to the thermal noise in terms of the dimen- 
sionless parameter 



a 



where (x^) 



ml 

2D 



1, andae(0,°o), (22) 



/eq 



and (x^)eq are the variance of x in the nonequilib- 
rium steady state (for /o ^ 0) and in equilibrium (for /o = 0) 
respectively. Using this parameter a and the above transfor- 
mations, the expression of Tj(A) given by ( [T3] l becomes 



Tj(A) = v/l+4aA(l-A). 
We also rewrite g{X) as 
2 



2tj(A) 



l + 77(A)-2aA l + Tj(A)+2aA^ 



(23) 



(24) 



where, the first factor in the above equation is due to the aver- 
aging over the initial conditions with respect to the the steady 
state distribution and the second factor is due to the integrating 
out of the final degrees of freedom. 

The probability density function of the work done is related 
to its characteristic function by the inverse Fourier transform 



1 

iTti 



8iX)e'f"^'^^dX, 



where 



/.,(A) = -[1-T7(A)]+Aw 



(25) 



(26) 



and we have set = 1 — this is equivalent to redefining the 
time in the unit of Ty, i.e., t/Tj^ t. The integration is done 
along the imaginary axis (vertical contour through the origin) 
in the complex-A plane. 



A. The large deviation function 

The large-T behavior of P{Wt) can be obtained from the 
saddle point approximation of the integral in Eq. p5] l. The 
saddle-point A* is obtained from the solution of the condition 

/;(A*)=Oas 



X*{w)^\ 



w 



1 + 1 
a 



(27) 



It follows from the above expression that X*{w) is a mono- 
tonically decreasing function of w (see Fig. [T]i and X*{w -> 
=Foo) A±, where 



1±^ 1 



(28) 



Therefore, A* e (A_,A+). The a dependence of A± are dis- 
played in Fig. [2] We note that r\ (A ) can be written in terms of 
A± as 



Tj(A) = v/4a(A+-A)(A-A_). 



(29) 



Clearly, Tj(A) has two branch points on the real-A line at A±, 
and Tj(A) is real and positive for A G (A_,A+). Therefore, 
/vi,(A) is also real on the real-A line in the interval (A_,A+). 
At A = A* we find 



T7(A*) = 



Vo:(l + «) 

^/w^ + a 



(30) 



Let us now look at the analyticity of ^(A ) given by Eq. ( 24 1 
for A e (A_, A+). It is clear that g{X) is real for A e (A_, ATJ. 
From Eq. ([28]l we find that 



2aA_ + l = Vl + a[\/l + a-\/a| >0. 



(31) 



Therefore, 1 + 77(A) +2aA > for A e (A_, A+) for afl a e 
(0,00). Again, from Eq. (28 1 we find that 



2aA+ - 1 = y/a{l + a) - (1 - a). 



(32) 
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FIG. 1. (Color online). Plots of the saddle-point X* as a function of 
w for a = 3/2 (red), 3 (blue) and 10 (orange). The horizontal dashed 
lines mark the positions of the pole Xq. The vertical dashed lines 
mark w* for respective values of a. 
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FIG. 2. (Color online). A+, and Xq are plotted (in red, orange, 
and blue colors, respectively) against a. 



The right side of the above equation is negative for a < 1/3. 
Thus, 1 — 2aX > for A < A+ and a < 1/3. Consequently, 
when a < 1/3, we have 1 + 77(A) — 2aA >OforA e (A_,A+). 
Therefore, g{X) does not have any singularities in the interval 
(A_,A+) as long as a < 1/3. Ignoring the subleading contri- 
bution g{X) in the saddle-point calculation gives 



WT] 



with 



hi{w) = MX*). 



(33) 



(34) 



By substituting the expression of A* in /,y(A*), after some 
algebra, we find that 



hi (w) 



- w - 



a 



(35) 



We now consider the case a> 1/3. In this case, it is useful 
to express g{X) in the form 



T7(A)[Tj(A) + 2aA-l] 1 



where 



a(l + a)A[l + 77(A) + 2aA] A-Aq ' 



Ao = 



l + a 



(36) 



(37) 



and A_ < < Ao < A+ (see Fig.[2]i. In order for Aq to be a zero 
of 1 + Tj (A) — 2aA, it has to satisfy the condition 2aAo — 1 = 
> 0. Using the above expression of Aq we get 



77(Ao) =2aAo- 1 = 



3a -1 
l + a ■ 



(38) 



Therefore, ^(A) possesses a simple pole at Aq when a > 1/3. 
As shown in Fig.[T[ for any given a, as w decreases from +°o 



to — oo, the saddle point A*(w) moves unidirectionally from 
A_ to A+ on the real A line. For sufficiently large w, we have 
A- < A* < Ao. In such a situation, the contour of integration 
can be deformed smoothly through the saddle point A* and 
therefore we still have P{Wz = wt) ^ e^''' . However, as one 
decreases w, at some particular value w = w*, the saddle-point 
hits the singularity — where w* is found by solving A* (w*) = 
Ao as 



a(a - 3) 
3a -1 



(39) 



For w < w* , we then have < Aq < A * . Therefore, while shift- 
ing the contour of integration in the complex-A plane, from 
its original path along the imaginary-A axis to the steepest de- 
scent path through the saddle -point A*, a contribution from 
the pole is picked up according to the residue theorem, which 
to the leading order is exp[T/„, (A())]. Now the second deriva- 
tive of fw{X) along the real-A axis at A* can be found to be 



2(w2 + a)3/2 



(40) 



The above expression is always positive, which means that 
fw{X) has a minimum at A* along real-A. That implies 
fw{X*) < /n.(Ao). Consequently, the leading saddle-point 
contribution exp[T/,v(A*)] is smaller than the leading pole 
contribution. Therefore, the leading behavior of P{Wx) is 
given by 



P{Wt: = wt) - e 



T/72(vi 



(41) 



with 



/i2W=/w(Ao). (42) 
By substituting Ao in fn.(Xo), after some algebra, we find that 

1 - a 2w 



/l2(w) 



l + a l + a' 



(43) 



FIG. 3. (Color online). P(Wt) against the scaled variable w = W-c/t 
for T = 100 and a = 3. The points (blue) are obtained from numerical 
simulation, and the solid line (red) plots the large deviation form 
given by Eq. |47j. The dashed (magenta) line plots the asymptotic 
form given by Eq. (jSTJ. 
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Re(A) 

FIG. 4. (Color online). The paths of the steepest descent for w > 
(red), w = (blue) and w <0 (orange) respectively. 



To summarize, the large deviation function, defined by 



h{w) = lim - \nP{Wr = wt), 

is given by the following. For a < 1/3: 

h{w)=hi{w), 

and for a > 1/3: 



h{w) 



hi{w) for w > w* 
/i2(w) for w <w* 



(44) to deform the original contour of integration in Eq. ( 25 1 into 



(45) 



(46) 



where hi{w) and h2{w) are given by Eq. (35i and Eq. (43 1, 
respectively, and w* is given by Eq. (39i. It is easy to show 
that/ii(w*) =/z2(w*) and/z'i(w*) 



the steepest descent path that passes through X* and the imag- 
inary part of fw{X) is constant along the new path. Setting 
Im[/w(A)] = constant = lm[ f„{X*)\ = gives the path of the 
steepest descent as (see Fig.|4|i 



w 
/a 



a 



a 



-4A2 



(48) 



where = Re(A) and A/ = lm(A) are, respectively, the real 
and imaginary parts of A, that is, X — Xr + iXj. The steepest 
descent path intersects the real-A axis at A* at an angle (p = 
n/2, which is evident from Eq. (40l as well as Eq. (48 1. 



B. Finite time corrections 



In the previous subsection we have obtained the leading be- 
havior of P(Wt), which has the large deviation form 



tA(vi 



(47) 



P{Wt = wt) - e 



with the large deviation function given by either Eq. ( [45[ ) or 
Eq. (46 1 depending on whether a < l/3ora> 1/3. In Fig. [3] 
we plot this form together with P{Wt:) that we have obtained 
from numerical simulation for T = 100. It is evident from the 



comparison that the large deviation form given by Eq. (47 1 is 



not adequate to explain simulation data (or experimental data) 
of finite time. In this subsection we obtain the sub-leading 
contributions to the large deviation form. 

We have found in the previous subsection that the saddle- 
point is located at A* e (A_ , A+) given by Eq. ( 27 1. We intend 



1. Case: a < 1/3 



We have found in Sec. Ill A that g{X) does not have any 
singularities in (A_,A+) for a < 1/3. Therefore, we can 
smoothly deform the contour along the path of the steepest 
descent through the saddle-point. Subsequently, following the 
the usual saddle-point approximation method, we write 



V2ng{X*)e^f-^^*^e'^ 



(49) 
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Using Eqs. ([30|, and with 

some algebra, we find that 



K{w) 



^ V2g{X*) 

l + (w + a)A*(w)-/ii(w) 



{/!2 (w) - /!i (w) } + (w - a) { A * (w) - Ao } 



(50) 



It should be noted that, although the last line in the above ex- 
pression is written in that particular fashion involving Ao and 
/!2(w), instead of [1 + (vv — a)X*{w) — it does not 
have any singularities for a < 1/3 as we have discussed in the 
previous subsection. Even for a = 1/3, the above expression 
diverges only in the limit of w ^ —°o. 

Finally, using Eq. ^9\ , and Eq. ( [50| , in ( p5| ), we find that 
(recall that <j) ^ 7t/2) 



TIT 



(51) 



It is seen from Fig. |5f a) that the above form is a perfect fit for 
the data obtained using numerical simulation for a = 1/3. 



2. Case: a > 1/3 

L et us consider the case a > 1/3. We have found in 
Sec. Ill A that g{X) for this case has a pole at Ao e (0,A+] 



and A* (w) < Ao for w>w* , whereas X*{w) > Aq for w < w* . 
Now if w ^ w*, then the asymptotic behavior of Eq. (25 i is 
obtained from the saddle-point approximation, which is given 
by Eq. (51 1. On the other hand, when w ^ w*, while deform- 



ing the contour of the integration in Eq. ( 25 1 to the steepest 
descent path through the saddle point, it goes around the pole 
at Ao in a clockwise direction. Therefore, according to the 
residue theorem, the integral pick up a contribution from the 
pole which is given by — ^-i exp[T/H.(Ao)], where 



: Um [(A 



'Ao)g(^)] 



(3a -1)^ 
8a2(l 



-a) 



(52) 



It is convenient to use Eq. ( 36 1 for g(X) to obtain the last ex- 
pression. The saddle-point contribution to the integral is same 
as given by Eq. (51 1. Therefore, combining the contributions 



from the saddle-point as well as from the pole, we can write 
that for a > 1 /3 when |w - w* | > 0, 



(53) 



Since h\{w*) = /i2(w*) and A*(vv* 
Eq. (50 1 that K(w) diverges as w — > 



Ao, it is clear from 
. Therefore, Eq. (53 1 



does not provide the correct description of the actual proba- 
bility density function as w approaches w* (from any side), 
which is also seen from Fig.|5]^b) and Fig.|5]^c). 

We carry out an asymptotic analysis of the integral along 
the path of the steepest descent, which is also valid when w 
is near w*, using the method of uniform asymptotic expan- 
sions |40|. The steps are outlined in Appendix [B| following 
which we get 



(»■) 



K(w) 



sgn(w* -H')g-i 

^/h2{w)-hl{w) 



r-1 



sgn(w* — w) 



erfc (v/t[/z2(h')-/ii(w)]) - {w* - w) 



(54) 



Again, this above asymptotic form matches with simulation 
results extremely well even for T = 10 as seen from Fig.|5|b) 
and Fig.|5|c). It follows from Appendix[B|that the above equa- 
tion reduces to Eq. (53 i for |w — w* | ^0. Equation ( [54) i is 
valid for any values of w including w = w* . 



IV. REMARKS ON STOCHASTIC INTEGRATION 

Finally, we make some remarks on the evaluation of the 
integrals of type Jq f{t)v{t)dt — that appears in Eqs. (|5]) and 
( pOj l — in numerical simulations or from experimental data. 
First the total time interval (0,t) is divided into small time 
intervals of Af duration such that 



Jo 



f{t)v{t)dt 



N-l 

L 

n=0 



r+Al 



f{t'W)dt' , 



(55) 



where t = nAt and NAt = T. Next the integral on the right side 
of the above equation needs to be evaluated for each interval 
(f ,? + Af). At first glance, it might look reasonable to assume 
that the integrals could be evaluated numerically in simula- 
tions or from experimental data by any one of the following 
approximation schemes: 



(1) 



(2) 



t+Ar 



fit'Ht')dt' 



''v{t)Ft^{tl 
1 



j'^"^ f{t')v{t')dt' « ^-[v{t)+v{t+At)]F^{t), 

rt+At 



(3) J f{t')v{t')dt' ^v{t + At)Fi^{t), 

where Ft^{t) = J,'^^ f{t')dt' , which is of 0{VAt). However, 
it turns out that the second one (a la Stratonovich) is the only 
correct scheme to follow. If one were to use the first scheme 
instead, the resulting distribution would shift to the left of the 
true one, whereas the use of the third scheme would result 
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FIG. 5. (Color online). /"(Wr) against the scaled variable w = IVt/t for T = 10 and (a) a = 1 /3, (b) a = 3, and (c) a = 20, respectively. 
The points (blue) are obtained from numerical simulation, and the dashed solid lines (magenta) plot the analytical asymptotic forms given by 
Eq. \5l\ for (a) and Eq. ^54\ for (b) and (c). In (b) and (c): the solid (red) lines plot the form given by Eq. ^53\ , and the dashed vertical lines 
mark the positions of w = w*. 
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1:^ 10-4 
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q;=18.66 
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FIG. 6. (Color online) The three sets [(1) red, (2) green and (3) blue] 
of points plot the probability density functions of the work fluctua- 
tions obtained from numerical simulations of the Langavin equation 
given by Eq. ^19\ , while using the three different numerical integra- 
tion schemes respectively, outlined in Sec. |IV| for the work integral 
given by Eq. ^20\ . The dashed line (magenta) plots the analytical 
asymptotic forms given by Eq. \54\ . It is evident that only the second 
stochastic integration scheme of Sec. |IV| yields the correct distribu- 
tion. 



in a shift to the right of the actual distribution (see Fig. |6|. 
While comparing our asymptotic form Eq. ( 54 1 with the ex- 
perimental distribution of Ref. [33], the work is evaluated us- 
ing the second scheme given above jJTl . The second scheme 
is also used while deriving the Fokker-Planck equation given 
by Eqs. (|7| and ([8]). On the other hand, the methods used here 
to obtain the solution of the Fokker-Planck equation does not 
require any such stochastic integration scheme 1 36 1 . We have 
used the Fokker-Planck equation to merely verify the solution. 
This also proves the correctness of the second scheme in the 
present context. 



V. SUMMARY 

In this paper, we have applied the recent formalism of 
Kundu et al. |36| to the case of a harmonically bound Brow- 
nian particle coupled to two heat baths at different tempera- 
tures. We have considered the fluctuations of the total amount 
of heat flow Q from one of the baths to the particle in a 
given time duration T. Its characteristic function for given 
initial and final phase-space configurations satisfies a non- 
trivial Fokker-Planck equation. We have obtained the largest 
eigenvalue of the Fokker-Planck operator as well as the corre- 
sponding eigenfunctions (left and right) exactly. Using those, 
and integrating over the final configurations and averaging 
over the initial configurations with respect to the nonequilib- 
rium steady-state, we have obtained the characteristic func- 
tion (e^^Q) w ^(A) exp[T/i(A)] for large T, where the both 
/i(A) and g{X) have been obtained exactly. A special case 
of this problem coiTesponds to the work done by an external 
stochastic force on a harmonic oscillator coupled to a ther- 
mal bath. This special case in fact models an recently studied 
experimental system of a stochastically driven atomic-force 
microscopy cantilever |33|. For this model we have analyzed 
the characteristic function and found the exact large devia- 
tion function as well as the complete asymptotic forms of the 
probability density function of the work. We have compared 
the analytical results with numerical simulations and found 
excellent agreements between the theory and simulation. In 
fact, the theoretical asymptotic forms of the probability den- 
sity function also match quite well with the experimentally 
obtained forms |37|. Finally, we believe that the results as 
well as the analytical methods of this paper would be useful 
for many other similar problems. 
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Appendix A: Fluctuations of the heat flow in a harmonic chain 

Here, for ease of reference, we outline the relevant re- 
sults, of Ref. |36|, for a harmonic chain consisting of 
particles connected at its two ends to white-noise Langevin 
reservoirs at different temperatures Tl and Tr respectively. 
The system, described by the variables = {xi,X2, ■ ■ ■ ,xn) 
and = (vi, V2, . . . , v^), evolves according to the following 
equations of motion: 



x = Afy = -4>x-yy + Tj(f), 



(Al) 



where M = diag(mi,m2, . . . ,mAr) is the mass matrix, 4> is 
the force matrix, y is the dissipation matrix with elements 
Yj j = 5/.y(5,',i7z, + 5,_Ar7s), and rj is the Gaussian noise vector 
with elements Tj, (f) — 5, j + 5,,AfTj^(r) whose correlators 
are given by Eqs. (|2])-(|4]i. The quantity of interest is the total 
amount of heat flowing from one of the reservoirs — say the 
left (L) — into the system in a given time duration T, given by 

Q = l' [ridt) - Yin (0] vi (?) dt. (A2) 

It was found in Ref. |!36| that the restricted characteris- 
tic function Z{X,U, tpo) ^ (e'^^ d[U - U{t)])jj^ for fixed 

initial and final configurations, Uq = {Xq,Vq ) and = 
{X^ ,V^), respectively, has the large-T asymptotic form 

Z{X,U,z\Uo) - x{Uo,X)^'{U,X) exp[TM(A)]. (A3) 

The cumulant generating function (A ) is given by 

1 



tl(A) = / dcoln 

47tJ- 



l+X{AI5-X)^i{co) 



(A4) 



where A)3 (jR/dR) - (YL/di) and 

^i{co)=4dLdRO)^Gl^{o))Gl^ico), 

with 



G^{co)= ^-co^M±icoY 



-1 



(A5) 



(A6) 



The above formula of /i(A) has been also generalized to 
the case of heat conduction across arbitrary harmonic net- 
works 1391 . 

The functions M^(f/,A) and x{Uq,X) have the following 
Gaussian forms: 



X(t/o,A) =exp 



exp[-if/^Li(A)f/] 



--U^L2[X)Uo 



where 



Lr{X)=H-'+H~,'Hl, 
L2{X) = -H-'Hl 

and Hi and H2 satisfy the relation 



(A7) 
(A8) 



(A9) 
(AlO) 



(All) 



The matrices Hi, H2, and H3 are, respectively, given as fol- 
lows: 



Hi(A) = i(/i+/[) with /i(A) = ^ 

2 7Z 



dco 



Ci,iF2F] +Ci,2F3F^ +C2,iF2F^ + C2,2F3F^ 



l+A(A/3-A)^i(a)) 
diil- 2icoYLG+i )F{F^ - 2ico{YL + XdL)dR G+f^F*F^ 



l+A(Aj3-A).5^i(a)) 



H2(A) = lim- / dcoe" 

1 T X(yt +Xdi) 

H3(A) = -(/3+/3 ) with /3(A) = ^^^^ ^/ dco 



FiF 



n J-^ l+A(Aj3-A)J^i((o) 
In the above expressions Fi, F2, and F3 are column matrices, given respectively by 

= ([G+4>] u , [G+4>] 1,2, [G+4>] 1,^, -/a)[G+M] la , -/a)[G+M] 1,2, . . 

P^I — ^ii' ■••'^7v,i' '(oG^^, icoG2 iT--,i(oG^ I ^, 

fJ' = (gIj^, G^j^, icoGlfj, '(^G^,N^ ■ ■ ■ ' ''^^n,n) ' 
and Cij is the {/, y'}-th element of the matrix 



/a)[G+M]i,^ , 



Cico) 



.2|^+ |2 



V 



--4XylCo'\GI^^ 
ur 

4A7LC(j2Gf J - 2/A coG 



1,N 



4XYLCO^G+fiij^ + 2iXcoGij^ 
^+4Xyr(o^\G+^\^ 



(A12) 
(A13) 
(A14) 

(A15) 
(A16) 
(A17) 



(A18) 
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1. The Af = 1 case 

Here we explicitly evaluate the above expressions for the 
special case of the single Brownian particle in a harmonic po- 
tential that is discussed in Sec. Ill] For this case, we have 



which gives 



G±(«) 



k— (0 m±iC0Y 



-1 



AdidRCO 



(A19) 



(A20) 



In this case Eq. \AA) can be evaluated explicitly and the 
final expression is given by Eq. ( 10 1. 



To obtain the explicit forms of of the functions given by 
Eqs. ( |A7| i and ( |A8| l, we first find that 



Fi=G^ 



k 

—iCOm 



and F2= = G 



(A21) 



Moreover, adding the elements of C{co) from Eq. ( A18 i and 
then using the identity —i[G^ — = 2y(o\G^ ^ we get 

Next, we carry out the integrations in Eqs. ( |A12[ )-( |AT41 ) to 
find that 



H2(A) 

^3(A) = 



dL^dR ( l/k 



y^{X) \ l/m J ' 
l£/i-(7/2)[Tj(A)-l] / 1 



77] (A) 

XiYL+Xd^) fk 

777(A) V ^ 



1 



(A22) 
(A23) 
(A24) 



It can be checked that the above matrices satisfy the relation 



given by Eq. (All 



Using the above matrices in Eqs. \A9\ and ( |A10| i, respec- 
tively, we obtain 



Li(A)=B+(A) 
L2(A)=B_(A) 



k 

m 

k 

m 



where 



7T?(A)±(7+2Afl'z.) 
2{dL + dR) 



(A25) 
(A26) 

(A27) 



Finally, using the above results in Eq. (jATji and Eq. ( A8 1 
we obtain Eq. (Ill and Eq. (12 1 respectively. 



Appendix B: Uniform asymptotic expansions for a saddle-point 
near a pole 

Following Ref. ll40l . we outline below the method of uni- 
form asymptotic expansion of a integral having a saddle-point 



near a pole. Let us denote 



(Bl) 



where C denotes the contour along the path of steepest de- 
scent. We set 



/„.(A)-.a,(a*) = -m2 



(B2) 



Since the imaginary part of /w(A) is constant (which happens 
to be zero in our particular case) on C, the above integral can 
be converted into an integral with respect to a real variable u, 

as 



/(T) 



q{u)i 



du, 



with 



dX 
du 



(B3) 



(B4) 



The pole Aq of ^(A) is then mapped to a pole of qiu) at m = ih 
with 



h = sgn(A* - Ao) - /.v(A*) 



(B5) 



Note that, /n(A) is minimum at A* along real-A. Therefore, 
/w('^) > /w('^*)i and hence h is real. We can write 



q{u) = 



c-l 
u — ih 



\j/{u), 



(B6) 



where c_i = lim„^,i[(M — //7)^(m)], and yf(u) has no pole 
at M = it or u — 0. To evaluate c_i we note that + 
= /n,(Ao) -/,v(A) and dX/du = -2m//' (A). Moreover, 

fwiX) = /vt.(Ao) + (A - Ao)/'(Ao) H as A Aq. Therefore, 

we get 



c_i = lim 



u + ib 



dX 
du 



: lim [(A-Ao)g(A)] 



(B7) 



where g^i is given by Eq. ( 52 1. 
Let us first look at the integral 



- du. 



/-oo u — ib 
It satisfies the differential equation 

j[{T)^b^JliT) = -ib^{^), 



(B8) 



(B9) 



with the boundary condition Ji {°o) = 0. It can be verified that 
the solution is given by 



7i(t) = sgn{b)i7:e'^''^eTfc{y/t\b\). 
Next, we look at the integral 

72(t) = / \j/{u)e' 



-^"^du. 



(BIO) 



(BID 
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Expanding yf{u) Taylor series about u = Q yf{u) = 
I^r=o (0)m"/"' ill the above integral, and then term by 
term integration of the series yields 



^(2m)(o) r(m+l/2) 



^0 (2m)! T'«+V2 ■ 
In terms of Ji.ii'^) we finally have 

/(T)=e^/"'(^*)[g_iyi(T)+72(T)]. 



(B13) 



In the following, we consider the large T limit and find the 
leading order contribution of /2(t), i.e.. 



72(T) = r(0)y(^+O(T-3/2). 



From Eq. ( B6 1 we get 

V^(0)-?(0) 
and then from Eq. (|B4|i 



8-1 
ib ' 



(B14) 



(B15) 



?(0)-g(A*) 



du 



(B16) 



Near A* we have 



*\2 



-u'^-f"{l*){X-X*) 



(B17) 



Let A - A* = re'^ and /"(A*) = |/"(A*)|e'^ Therefore, 



-|/"(A*)|e'('^+2^)- 



(B12) Since, m is real, we have 9 + 2(j) = n, and 



±-^|/"(A*)|i/2- 



Thus, 

du 
dX 



1— j-A,* 



£/m dr 
dr dX 



V2 



1/2 .-'-0 



Substitution of this into Eq. ( |B16[ ) yields 



(B18) 



(B19) 



(B20) 



(B21) 



In our case, from Eq. (40i we get 0=0, and therefore = 
n/2. 

If the saddle- point and the pole are far apart, then using 
Ji{x) ^ {i / b) yjjnjx) for large b and T, it is immediately 
checked that /(t) reduces to the usual saddle-point approx- 
imation given by Eq. (49 1. More generally, /(t) given by 
Eq. (B13 I, with 71,2(7) given by Eqs. ( |B10| l and ( |B12| ) respec- 
tively, IS valid for all b, including the limit — > 0. For large T, 
using Eqs. ( |B14| l, ( |B15| l and ( |B21| i we get 



-/2(T) 



1% 



V2g{X*) 



(B22) 
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